Building intuition into popgen fundamentals and inference

Workshop on population and speciation genomics

Many problems in population genetics cannot be solved by a mathematician, no matter how gifted. [It] is already clear that computer methods are very powerful. This is good. It […] permits people with limited mathematical knowledge to work on important problems […].

James F. Crow – interview

Why use simulations?

1. Developing intuition into statistics
2. Estimating model parameters (i.e. ABC)
3. Ground truth for method development

Developing intuition into statistics

Image from Peter (2016)

Developing intuition into statistics

Image from Lawson et al. (2018)

Estimating model parameters (i.e. ABC)

Image from Wikipedia on ABC

Ground truth for method development

Image from Schiffels and Durbin (2014)

Simulation software

The most famous and widely used are SLiM and msprime.

They are very powerful but both require:

• quite a bit of programming knowledge,
• a lot of code for non-trivial simulations (🐛🪲🐜).

Our exercises will focus on the slendr simulation toolkit for population genetics in R.

But first let’s look at SLiM and msprime at least a little bit…

What is SLiM?

• A forward-time simulator

• Has its own programming language!

• Massive library of functions for:
  • demographic events
  • various mating systems
  • natural selection

• More than 700 pages long manual!

Figure modified from Alexei Drummond

Simple neutral simulation in SLiM

initialize() {
    // create a neutral mutation type
    initializeMutationType("m1", 0.5, "f", 0.0);

    // initialize 1Mb segment
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, 999999);

    // set mutation rate and recombination rate of the segment
    initializeMutationRate(1e-8);
    initializeRecombinationRate(1e-8);
}

// create an ancestral population p1 of 10000 diploid individuals
1 early() { sim.addSubpop("p1", 10000); }

// in generation 1000, create two daughter populations p2 and p3
1000 early() {
    sim.addSubpopSplit("p2", 5000, p1);
    sim.addSubpopSplit("p3", 1000, p1);
}

// in generation 10000, stop the simulation and save 100 individuals
// from p2 and p3 to a VCF file
10000 late() {
    p2_subset = sample(p2.individuals, 100);
    p3_subset = sample(p3.individuals, 100);
    c(p2_subset, p3_subset).genomes.outputVCF("/tmp/slim_output.vcf.gz");

    sim.simulationFinished();
}

What is msprime?

• A Python module for writing coalescent simulations

• Extremely fast (genome-scale, population-scale data)

• You must know Python fairly well to build complex models

Figure modified from Alexei Drummond

Simple simulation using msprime

import msprime

demography = msprime.Demography()
demography.add_population(name="A", initial_size=10_000)
demography.add_population(name="B", initial_size=5_000)
demography.add_population(name="C", initial_size=1_000)
demography.add_population_split(time=1000, derived=["A", "B"], ancestral="C")

ts = msprime.sim_ancestry(
  sequence_length=10e6,
  recombination_rate=1e-8,
  samples={"A": 100, "B": 100},
  demography=demography
)

And yet I built something else 🙄

www.slendr.net

… but why?

First motivation: spatial simulations!

A somewhat broader motivation

• Most researchers are not expert programmers

• All but the most trivial simulations require lots of code

• Yet, 90% ^{[citation needed]} of simulations are basically the same!

  • create populations (splits and \(N_e\) changes)

  • specify if/how they should mix (rates and times)

  • … all this means duplication across projects

• Computing statistics presents even more hurdles

slendr makes this trivial, even for complex models

slendr crash course

Typical steps of a slendr R workflow

We will always start with this:

library(slendr); init_env()

The interface to all required Python modules has been activated.

Followed by some combination of the following:

1. creating populations
2. programming \(N_e\) size changes
3. encoding gene-flow events
4. simulating genomic data
5. computing popgen statistics

Creating a population

At minimum, we need its name, size and “time of appearance”:

pop1 <- population("pop1", N = 1000, time = 1)

. . .

This creates a normal R object. Typing it out gives a summary:

pop1

slendr 'population' object 
-------------------------- 
name: pop1 
non-spatial population
stays until the end of the simulation

population history overview:
  - time 1: created as an ancestral population (N = 1000)

Note: Because slendr uses either msprime or SLiM internally for simulation of genomic data, all individuals are assumed to be diploid.

Programming population splits

Splits are defined by providing a parent = <pop> argument:

pop2 <- population("pop2", N = 100, time = 50, parent = pop1)

Note: Here pop1 is an R object created above, not a string "pop1"!

. . .

The split is reported in the “historical summary”:

pop2

slendr 'population' object 
-------------------------- 
name: pop2 
non-spatial population
stays until the end of the simulation

population history overview:
  - time 50: split from pop1 (N = 100)

Scheduling resize events

• Step size decrease:

pop1 <- population("pop1", N = 1000, time = 1)
pop1_step <- resize(pop1, N = 100, time = 500, how = "step")

• Exponential increase:

pop2 <- population("pop2", N = 100, time = 50, parent = pop1)
pop2_exp <- resize(pop2, N = 10000, time = 500, end = 2000, how = "exponential")

Tidyverse-style pipe %>% interface

• Step size decrease:

pop1 <-
  population("pop1", N = 1000, time = 1) %>%
  resize(N = 100, time = 500, how = "step")

• Exponential increase:

pop2 <-
  population("pop2", N = 1000, time = 1) %>%
  resize(N = 10000, time = 500, end = 2000, how = "exponential")

Note: You can read (and understand) a() %>% b() %>% c() as “take the result of the function a, pipe it into function b, and then pipe that to function c”.

A more complex model

pop1 <- population("pop1", N = 1000, time = 1)

pop2 <-
  population("pop2", N = 1000, time = 300, parent = pop1) %>%
  resize(N = 100, how = "step", time = 1000)

pop3 <-
  population("pop3", N = 1000, time = 400, parent = pop2) %>%
  resize(N = 2500, how = "step", time = 800)

pop4 <-
  population("pop4", N = 1500, time = 500, parent = pop3) %>%
  resize(N = 700, how = "exponential", time = 1200, end = 2000)

pop5 <-
  population("pop5", N = 100, time = 600, parent = pop4) %>%
  resize(N = 50, how = "step", time = 900) %>%
  resize(N = 1000, how = "exponential", time = 1600, end = 2200)

Again, each object carries its history!

For instance, this is the summary you will get from the last population from the previous code chunk:

pop5

slendr 'population' object 
-------------------------- 
name: pop5 
non-spatial population
stays until the end of the simulation

population history overview:
  - time 600: split from pop4 (N = 100)
  - time 900: resize from 100 to 50 individuals
  - time 1600-2200: exponential resize from 50 to 1000 individuals

. . .

This way, you can build up complex models step by step, checking everything by interacting with the R console.

Gene flow / admixture

We can schedule gene flow from pop1 into pop2 with:

gf <- gene_flow(from = pop1, to = pop2, start = 2000, end = 2200, rate = 0.13)

Note: Here rate = 0.13 means 13% migrants over the given time window will come from “pop1” into “pop2”.

. . .

Multiple gene-flow events can be gathered in a list:

gf <- list(
  gene_flow(from = pop1, to = pop2, start = 500, end = 600, rate = 0.13),
  gene_flow(from = ..., to = ..., start = ..., end = ..., rate = ...),
  ...
)

Model compilation

model <- compile_model(
  populations = list(pop1, pop2, pop3, pop4, pop5),
  generation_time = 1,
  simulation_length = 3000,
  direction = "forward"
)

This takes a list of model components, performs some consistency checks, and again returns a single R object.

Model compilation

model <- compile_model(
  populations = list(pop1, pop2, pop3, pop4, pop5),
  gene_flow = gf,      # <----- in case our model includes gene flow(s)
  generation_time = 1,
  simulation_length = 3000,
  direction = "forward"
)

Gene flow(s) are included via the gene_flow argument.

Model summary

Typing the compiled model into R prints a brief summary:

model

slendr 'model' object 
--------------------- 
populations: pop1, pop2, pop3, pop4, pop5 
geneflow events: 1 
generation time: 1 
time direction: forward 
time units: generations 
total running length: 3000 time units
model type: non-spatial

configuration files in: /private/var/folders/70/b_q2zdh116b9pfg29p03sx600000gn/T/RtmpmR3uVo/file1843e774e3aa7_slendr_model 

This can be useful as a quick overview of the model we are working with. However, a better way to check a model is…

Model visualization

plot_model(model)

A note on units (and direction) of time

“Forward time units”

pop1 <- population("pop1", N = 1000, time = 1)

pop2 <-
  population("pop2", N = 1000, time = 300, parent = pop1) %>%
  resize(N = 100, how = "step", time = 1000)

pop3 <-
  population("pop3", N = 1000, time = 400, parent = pop2) %>%
  resize(N = 2500, how = "step", time = 800)

pop4 <-
  population("pop4", N = 1500, time = 500, parent = pop3) %>%
  resize(N = 700, how = "exponential", time = 1200, end = 2000)

pop5 <-
  population("pop5", N = 100, time = 600, parent = pop4) %>%
  resize(N = 50, how = "step", time = 900) %>%
  resize(N = 1000, how = "exponential", time = 1600, end = 2200)

model <- compile_model(
  populations = list(pop1, pop2, pop3, pop4, pop5),
  generation_time = 1,
  simulation_length = 3000,
  direction = "forward"
)

Note: We start pop1 in generation 1, and each later event has an increasing time value.

“Forward time units”

plot_model(model, proportions = TRUE)

Note: We start pop1 in generation 1, and each later event has an increasing time value.

“Backward time units”

pop1 <- population("pop1", N = 1000, time = 30000)

pop2 <-
  population("pop2", N = 1000, time = 27000, parent = pop1) %>%
  resize(N = 100, how = "step", time = 20000)

pop3 <-
  population("pop3", N = 1000, time = 26000, parent = pop2) %>%
  resize(N = 2500, how = "step", time = 22000)

pop4 <-
  population("pop4", N = 1500, time = 25000, parent = pop3) %>%
  resize(N = 700, how = "exponential", time = 18000, end = 10000)

pop5 <-
  population("pop5", N = 100, time = 24000, parent = pop4) %>%
  resize(N = 50, how = "step", time = 21000) %>%
  resize(N = 1000, how = "exponential", time = 14000, end = 8000)

model <- compile_model(
  populations = list(pop1, pop2, pop3, pop4, pop5),
  generation_time = 10 # (10 time units for each generation)
)

Note: Same model as before, except now expressed in units of “years before present”.

“Backward time units”

plot_model(model, proportions = TRUE)

Note: Same model as before, except now expressed in units of “years before present”.

So we built a model…

… but how do we simulate data from it?

Built-in simulation “engines”

slendr has two simulation “engine scripts” built-in:

• SLiM engine (slendr source)
• msprime engine (slendr source)

. . .

They are designed to understand slendr models, which means that you can simulate data just with this command:

ts <- msprime(model, sequence_length = 100e6, recombination_rate = 1e-8)

. . .

No need to write any msprime or SLiM code!

The result of a simulation is a tree sequence (ts)

What is tree sequence?

• a record of full genetic ancestry of a set of samples
• an encoding of DNA sequence carried by those samples
• an efficient analysis framework

Why tree sequence?

Why not VCF or a normal genotype table?

What we usually have

What we usually want

An understanding of our samples’ evolutionary history:

This is exactly what a tree sequence is!

The magic of tree sequences

They allow us to compute statistics without genotypes!

There is a “duality” between mutations and branch lengths.

See an amazing paper by Ralph et al. (2020) for more detail.

What if we need mutations though?

Coalescent and mutation processes can be decoupled!

This means we can add mutations to ts after the simulation using ts_mutate().

Let’s go back to our example model…

… simulate a tree sequence…

In our script we’ll have something like this:

library(slendr)
init_env()

<... population() definitions ...>

model <- compile_model(...)
  
ts <-
  msprime(model, sequence_length = 50e6, recombination_rate = 1e-8)

… and overlay mutations on it

In our script we’ll have something like this:

library(slendr)
init_env()

<... population() definitions ...>

model <- compile_model(...)
  
ts <-
  msprime(model, sequence_length = 50e6, recombination_rate = 1e-8) %>%
  ts_mutate(mutation_rate = 1e-8)

In some exercises, mutations won’t be necessary. Where we will need them, you can use ts_mutate() like this.

So we can simulate data

How do we work with this ts thing?

Standard genotype formats

If a tree sequence doesn’t cut it, you can always:

• export genotypes to a VCF file:

ts_vcf(ts, path = "path/to/a/file.vcf.gz")

• export genotypes in an EIGENSTRAT format:

ts_eigenstrat(ts, prefix = "path/to/eigenstrat/prefix")

• access genotypes in a data frame:

ts_genotypes(ts)

   pos pop1_1_chr1 pop1_1_chr2 pop1_2_chr1 pop1_2_chr2 pop1_3_chr1 pop1_3_chr2
1  972           0           0           0           0           0           0
2 1051           0           0           0           0           0           0

slendr’s R interface to tskit statistics

Allele-frequecy spectrum, diversity \(\pi\), \(F_{ST}\), Tajima’s D, etc.

Find help at slendr.net/reference or in R under ?ts_fst etc.

General patterns of tskit computation

1. For a function ts_stat() and a slendr tree-sequence object ts, we compute a statistic (getting a data frame) like this:

df <- ts_stat(ts, sample_sets = <vector of individual names>)

2. If ts_stat() works between/across groups of individuals:

df <- ts_stat(ts, sample_sets = <list of vectors of individual names>)

3. For something like the \(f_3\) statistic, the arguments must be more precisely specified:

df <- ts_f3(ts, A = <ind. names>, B = <ind. names>, C = <ind. names>)

Specifying individuals to compute on

We can get samples recorded in ts with ts_samples():

ts_samples(ts) %>% head(2)

# A tibble: 2 × 3
  name    time pop  
  <chr>  <dbl> <chr>
1 pop1_1     0 pop1 
2 pop1_2     0 pop1 

. . .

A shortcut ts_names() can also be useful:

ts_names(ts) %>% head(5)

[1] "pop1_1" "pop1_2" "pop1_3" "pop1_4" "pop1_5"

. . .

We can also get a per-population list of individuals:

ts_names(ts, split = "pop")

$pop1
[1] "pop1_479" "pop1_721" "pop1_631" "pop1_432" "pop1_12" 

All slendr statistics accept this as their second (and further) arguments.

This is modelled after the sample_sets= argument of the respective tskit Python methods (except that you use “symbolic names” of individuals, not tree-sequence node numbers).

Example: nucleotide diversity

Get list of individuals in each population:

samples <- ts_names(ts, split = "pop")

names(samples)

[1] "pop1" "pop2"

samples$pop1 %>% head(3)

[1] "pop1_1" "pop1_2" "pop1_3"

samples$pop2 %>% head(3)

[1] "pop2_1" "pop2_2" "pop2_3"

 

Compute nucleotide diversity:

ts_diversity(ts, sample_sets = samples)

# A tibble: 2 × 2
  set   diversity
  <chr>     <dbl>
1 pop1  0.000401 
2 pop2  0.0000827

Example: allele frequency spectrum

Get names of all individuals:

samples <- ts_names(ts)[1:5]
samples

[1] "pop_1" "pop_2" "pop_3" "pop_4" "pop_5"

Compute the AFS:

afs <- ts_afs(ts, sample_sets = list(samples))

afs[-1]

 [1] 3917 2151 1432  941  740  624  607  587  416  385

 

plot(afs[-1], type = "b",
     xlab = "allele count bin",
     ylab = "frequency")

Working with multi-population statistics

For multipopulation statistics, you need a (named) list of samples recorded in each population.

ts_names() has a split = "pop" option just for this:

samples <- ts_names(ts, split = "pop")

$<NA>
NULL

$<NA>
NULL

You can use this in place of sample_names in code like:

ts_diversity(ts, sample_sets = samples)

More information

• slendr paper is now in PCI EvolBiol

• documentation, tutorials is here

• GitHub repo (bug reports!) is here

---

• check out my new demografr inference package

---

• contact details at bodkan.net

Let’s get started

Resources and exercise materials

These slides and other materials are (and always will be) at:

github.com/bodkan/cesky-krumlov-2025

On the GitHub page, you will find links to:

• these slides (for reference during exercise)
• a version of these slides as a single-page document
• exercise materials (with solutions built-in)

Open the GitHub link so you have everything at hand later.